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Abstract. In this paper we consider denoising and 
inpainting problems for higher dimensional combined 
cyclic and linear space valued data. These kind of 
data appear when dealing with nonlinear color spaces 
such as HSV, and they can be obtained by changing 
the space domain of, e.g., an optical flow field to po¬ 
lar coordinates. For such nonlinear data spaces, we 
develop algorithms for the solution of the correspond¬ 
ing second order total variation (TV) type problems 
for denoising, inpainting as well as the combination of 
both. We provide a convergence analysis and we apply 
the algorithms to concrete problems. 

Keywords. Higher order total variation minimiza¬ 
tion, vector-valued TV, cyclic data, combined denois¬ 
ing and inpainting, cyclic proximal point algorithm. 

1 Introduction 

One of the most well known methods for edge¬ 
preserving image denoising is the variational ap¬ 
proach minimizing the Rudin-Osher-Fatemi (ROF) 
functional [70]. In its basic form, it deals with scalar 
data. Related variational approaches for vector space 
valued data have gained a lot of interest in the liter¬ 
ature and are still topic of ongoing research; we ex- 
emplarily refer to [12, 42, 61, 67] and the references 
therein. In this paper, we consider TV-type function¬ 
als incorporating first and second order differences for 
the nonlinear data spaces which combine vector space 
valued data —in the following called linear space data 
to avoid confusion— and vectors of cyclic data. In 
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these spaces, we deal with denoising and inpainting 
problems as well as simultaneous inpainting and de¬ 
noising problems. 

Image inpainting is a problem arising in many ap¬ 
plications in image processing, image analysis and re¬ 
lated fields. Examples are restoring scratches in pho¬ 
tographs, removal of superimposed objects, dealing 
with an area removed by a user, digital zooming as 
well as edge decoding. Principally, any missing data 
situation —whatever the reason might be— results in 
an inpainting problem. This is not restricted to 2D im¬ 
ages. Further examples are defects in audio and video 
recordings, or in seismic data processing. In this re¬ 
spect, also interpolation, approximation, and extrapo¬ 
lation problems may be viewed as inpainting problems. 
We recommend the survey [45] and Chapter 6 of [27] 
as well as [16, 18] for an overview on inpainting and for 
further applications. There are various conceptionally 
different approaches to inpainting, cf. [27, 45] and [18] 
which also includes some comparison. Among these 
are methods based on linear transforms from harmonic 
analysis such as curvelets and shearlets which are com¬ 
bined with a sparsity approach based on minimiza¬ 
tion on the corresponding coefficients [17, 32, 33, 51]. 
Other approaches are based on (often nonlinear) PDE 
and variational models, cf., e.g., [13, 19, 25, 26, 34, 
57-60, 74, 78]. 

In general, exemplar-based and sparsity-based meth¬ 
ods perform better for filling large texture areas, 
whereas diffusion-based and variational techniques 
yield better results for natural images. Among the 
variational techniques applied, total variation (TV) 
minimization is one of the prominent models. The 
minimizer of the corresponding TV functional yields 
the inpainted image. TV inpainting works well for 
elongated inpainting areas but has problems with 
larger gap connections. In such situations higher or¬ 
der, in particular curvature based, schemes perform 
better. 
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The first TV regularized model was proposed in [70] 
for denoising. It was first applied to missing data situ¬ 
ations/inpainting in [4, 24]. Further references for TV 
based image inpainting are [15, 28, 74]. In contrast to 
classical methods, the results are typically not over¬ 
smoothed; however, it is well known that these min- 
imizers very often show ‘staircasing’ effects, i.e. the 
result is often piecewise constant, although the under¬ 
lying signal varies smoothly in the corresponding re¬ 
gions. In order to avoid staircasing, higher order and, 
in particular, second order differences and derivatives 
(in a continuous domain setting), are often employed. 
References are the pioneering work [20] as well as [14, 
21, 23, 30, 31, 48, 53, 55, 56, 71-73]. TV functionals 
for linear space valued data were considered in [12] in 
the context of linear color spaces. The papers [63, 64] 
deal with denoising and inpainting in the RGB color 
space using linear combinations of first and second or¬ 
der terms. A total generalized variational model can 
be found, e.g., in [61]. In contrast to applying pure 
second order terms or linear combinations of first and 
second order terms, total generalized variational ap¬ 
proaches try to find some optimal balancing between 
the first and second derivatives. The advantage is that 
the related schemes better preserve the edge struc¬ 
tures. A detailed description may be found in [14]. 
The authors of [61] obtain a model for denoising lin¬ 
ear space valued color data. Second order total gen¬ 
eralized variation was generalized for tensor fields in 
[80]. 

However, in many applications, data having values 
in nonlinear spaces appear. Examples are diffusion 
tensor images [5, 65], color images based on non-fiat 
color models [22, 50, 52, 81] or motion group-valued 
data [69, 79]. Due to its importance, processing such 
manifold valued data has gained a lot of interest in re¬ 
cent years. To mention only some examples, wavelet- 
type multi scale transforms for manifold data have 
been considered in [44, 79, 83]. Statistical issues on 
Riemannian manifolds are the topic of [10, 11, 37, 
62, 66] and circular data are, in particular, consid¬ 
ered in [36, 49]. Furthermore, manifold-valued partial 
differential equations are studied in [29, 43, 77]. 

For TV functionals for manifold-valued data, an 
analysis from a theoretical viewpoint has been car¬ 
ried out in [39, 40]. These papers extend previous 
work [38] on S^-valued functions where, in particular, 
the existence of minimizers of certain TV-type ener¬ 
gies is shown. An algorithm for TV minimization on 
Riemannian manifolds was proposed in [54]. This ap¬ 
proach uses a reformulation as a multi label optimiza¬ 
tion problem with an infinite number of labels and a 
subsequent convex relaxation. An approach for lin¬ 
ear spaces using a relaxation of the label optimization 
problem as well was presented in [41]. First order TV 
minimization for -valued data has been considered 
in [75, 76]. In particular, these authors consider in¬ 


painting for -valued data in a first order TV setup. 
We proposed a different approach to first order TV 
minimization for manifold-valued data via cyclic and 
parallel proximal point algorithms in [85]. We estab¬ 
lished a second order setup for denoising -valued 
data based on cyclic proximal point algorithms in [6]. 
Inpainting for -valued data was considered in the 
authors’ conference proceeding [7]. 

Data consisting of combined cyclic and linear space 
components appear in various contexts. For example, 
such data appear when dealing with nonlinear color 
spaces such as HSV, HSL, HSI or HCL. Such data 
also appear in the context of optical flows. When con¬ 
sidering the flow vectors between consecutive images 
in polar coordinates which means separating magni¬ 
tude and direction, the resulting data takes its values 
in M X S^. In this context, this approach is natural 
and interesting and seems promising to improve the 
results obtained with the usual M^-valued approach. 
In principle, whenever data is given as vectors of polar 
coordinates, we are in the combined cyclic and real¬ 
valued setup of data in x considered in this 

paper. 

Contributions. In this paper, we consider inpaint¬ 
ing, denoising as well as combined inpainting and de¬ 
noising problems for combined cyclic and linear space 
valued images in x based on a second or¬ 

der TV-type formulation. We consider two variational 
models: the first model deals with the noise free situ¬ 
ation whereas the second one also considers the noisy 
case combining denoising and inpainting (including 
pure denoising by specifying the inpainting area as the 
empty set). In our nonlinear setting, these higher or¬ 
der approaches avoid unwanted staircasing effects as 
well. For combined cyclic and linear space data, we 
derive solvers for these variational problems based on 
cyclic proximal point algorithms. We provide a con¬ 
vergence analysis for both the noisy and the noise free 
model based algorithms developed in this paper. For 
both algorithms, we show the convergence to a mini- 
mizer under certain restrictions which are typical when 
dealing with nonlinear data. In particular, we assume 
the data to be dense enough meaning that they are lo¬ 
cally (and not necessarily globally) nearby. We apply 
our algorithms to denoising, inpainting and combined 
denoising and inpainting in the nonlinear HSV color 
space. Furthermore, we apply our algorithms for de¬ 
noising frames in volumetric phase-valued data - in 
our case, frames of a 2D film. Our approach is based 
on utilizing the neighboring k frames to incorporate 
the temporal neighborhood. The idea generalizes to 
arbitrary data spaces and volumes consisting of layers 
of 2D data. 

The novelties of the present work in relation to the 
authors previous work [6, 7] are as follows, (i) In con- 
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trast to [6, 7] we consider the more general, practically 
relevant, data spaces (S^)^ x here. In contrast 
to general manifolds, these product spaces still bear 
enough structure relevant for our purposes. We point 
out that both, the algorithmic and analysis part, are 
more involved than only component-wise considering 
the -valued or real-valued situation, (ii) Concerning 
the algorithmic part, we compute proximal mappings 
of constrained problems arising in inpainting situa¬ 
tions in this work. This is even new for data. In [7], 
we used a less natural projection approach generaliz¬ 
ing [6]. (iii) Concerning the analysis, we here include 
an inpainting setup. The conference proceeding [7] 
does not contain an analysis and [6] considers func¬ 
tionals for denoising valued data. A more detailed 
discussion may be found at the end of the paper. 

Outline of the paper. In Section 2 we introduce 
the variational models we consider for inpainting and 
denoising of combined cyclic and linear space data in 
this paper. We start with vector space data in Subsec¬ 
tion 2.1; then we define absolute differences for com¬ 
bined cyclic and vector space data in Subsection 2.2 
which allow us to derive the corresponding variational 
models for combined cyclic and vector space data. 
This is done in Subsection 2.3 for both inpainting 
noise free combined cyclic and vector space data as 
well as inpainting and denoising combined cyclic and 
vector space data. In Section 3 we develop algorithms 
for minimizing the variational models introduced pre¬ 
viously. These algorithms base on the cyclic proxi¬ 
mal point algorithm we present in Subsection 3.1. We 
present explicit formulas for the proximal mappings 
needed for inpainting in Subsection 3.2. Using these 
explicit representations, we derive a cyclic proximal 
point algorithm for inpainting both the noisy and noise 
free combined cyclic and vector space data in Subsec¬ 
tion 3.3. The convergence analysis of both algorithms 
is the topic of Section 3.4. Finally, in Section 4 we 
apply the derived algorithms to various concrete situ¬ 
ations. We consider denoising data living in the non¬ 
linear HSV color space in Subsection 4.1. Then we 
consider inpainting for noise-free as well as noisy data 
in such color spaces in Subsection 4.2. Finally, we ap¬ 
ply our algorithms for denoising frames of a -valued 
2D film in Subsection 4.3. 

2 Second order variational models for 
inpainting and denoising combined 
cyclic and linear space data 

In this section we derive models for denoising, inpaint¬ 
ing as well as simultaneous inpainting and denoising 
data having cyclic and linear space components. In 
Subsection 2.1, we first concentrate on introducing the 
considered models based on first and second order ab¬ 


solute finite differences restricting to the linear space 
setting. In Subsection 2.2 we obtain suitable defini¬ 
tions for absolute differences for combined cyclic and 
vector space data. In Subsection 2.3, we use these 
definitions to obtain inpainting and simultaneous in¬ 
painting and denoising models for combined cyclic and 
linear data. In particular, denoising is covered by con¬ 
sidering the empty set as inpainting region. 

Let us fix some notations. We denote by 
^0 := X images 

of size N hy which can also be seen as func¬ 
tions on the image domain Uq. For any subset 

U C Uo of pixel indices := Uo\U denotes the com¬ 
plement of Q in Uq. Each entry fij G T is called a 
pixel value, where T is some data space. Our main fo¬ 
cus is the (m + n)-dimensional space T := xW^. 

On a data space dx{x^y)^Di{x^y)^D 2 {x^y^ z)^ and 
Di^i{w,x^y, z) denote the distance on T, the abso¬ 
lute first, absolute second, and absolute second order 
mixed differences. If the space is clear from the con¬ 
text or the setting holds for all spaces, we omit the 
T. The elements x = (xi,...,Xn)^ G are col¬ 
umn vectors, where denotes the transposition. For 
two vectors x, ^ G we denote by (x, y) := x^y the 
inner product. By {x) 27 r ^ [~'7r,7r) we denote the ele¬ 
mentwise modulo operation, i.e. the unique value such 
that X = 27 rk + (x)27r, k G Finally, for matrices 
X = G we denote the columns by 


2.1 Inpainting and denoising vector space 
data 

The Rudin-Osher-Fatemi (ROF) functional [70] 

+ay2\\^^i,j\\^ a>0, 

hi hi 

is one of the most well known and most popular func¬ 
tionals in variational image processing. In its penal¬ 
ized form, it consists of two terms: the first term mea¬ 
sures the distance to the data / the second term is a 
TV regularizer where V denotes the discrete gradient 
operator, usually implemented as first order forward 
differences in vertical and horizontal directions. Both 
an isotropic version using the euclidean length or 2- 
norm and an anisotropic version employing the sum 
of absolute values or 1-norm for the second term are 
widely used. In this work we will restrict ourselves to 
the anisotropic version. In this form, the ROF Model 
is typically used for denoising purposes. To avoid the 
appearing staircasing effect, often higher order and, 
in particular, second order differences (respectively, 
derivatives, in a continuous domain setting) are em¬ 
ployed [14, 20, 21, 23, 30, 31, 48, 53, 55, 56, 71-73]. 
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Denoising. For pure denoising we consider the dis¬ 
crete second order TV-type functional 

Jix)=F{x-,f) + aTV^{x) 

+ /3TV2(:r) + 7TVi,i(x), 

where x, / are images defined on the image domainf^o 
denotes the image domain. The data values fij it¬ 
self live in a certain data space. Then the data term 
F{x; f) for given data / reads 

N,M 

Fix-J) = 2 V ( 2 ) 

^i=i 

where d is a distance on the data space. For data liv¬ 
ing in a vector space, d{fij^Xij) = \\fij — Xij\\ is an 
appropriate choice. In the pure linear space data situ¬ 
ation, the above quadratic data term (2) corresponds 
to a Gaussian noise model. For other types of noise, 
different data terms are more appropriate; e.g., for 
Laplacian noise the term F{x; f) = Ylfj=i 
is appropriate; cf. also [63]. 

The first order difference component (aTVi(x) is 
given by 

N-1,M 

aTVi(a;) = Q;i Di(xij,Xi+ij) 

N,M-1 

i,j=l 

+ ^ V Di{xij,Xi+ij+i) 

^ i,j=l 

^ i,j = l 

Again, for a vector space any norm of the ordinary ab¬ 
solute first order difference = \\xi^j — 

is an appropriate choice. The first order TV 
term incorporates horizontal, vertical and both diag¬ 
onal differences. The diagonals are incorporated to 
reduce unwanted anisotropy effects and are scaled by 
^ to take the length of the diagonal on the pixel grid, 
i.e. the distance of two pixels, into account. We note 
that J'{x) = F{x] f) + ((ai, 0 ^ 2 , 0, 0) TVi(x) is just the 
vector version of the anisotropic discrete ROF func¬ 
tional above. Using the notation 

D 2 {x,y,z) = \\x-2y + z\\ (4) 

and 

Di,i{x,y,u,v) = \\x - y - u + v\\, 

for a norm of the standard second order differences for 
vector space data, the second order difference compo¬ 
nent, consisting of a horizontal and vertical component 


f3Ty2{x) as well as a diagonal component 7 TViq(x), 
is given by 

N-1,M 

/5TV2(^) — ^ ^ -^2 (^i—l,jf5 5 

(5) 

N,M-1 ^ ’ 

+/32 -P 2 

7TVi,i(a:;)=7 Di^i{xij,Xi+ij,Xij+i,Xi+ij+i) 

(6) 

We note that similar to the diagonal differences in the 
TVi part the diagonal component 7 TVip(x), which 
is based on the mixed second order difference 
reduces unwanted anisotropy effects for the second or¬ 
der part. Actually, Di^i{xij,Xi^ij,Xij^i,Xi^ij^i) 
may be interpreted as a diagonal difference: av¬ 
eraging Xij^i and Xi^ij yields an estimate for 
the non-grid value m = “^^^ 2 + 1 / 27 + 1 / 2 ” • Then 

-^2 5'^5 ^i+l,j + l) — -^l,l(^i,jf 5 

X 2 +i,j, X 2 ,j+i, X 2 + 17 + 1 ) plays the role of a second order 
diagonal difference. Interchanging the roles of ^r^j+i, 
x^+ij and Xij,Xi+ij+i yields the other diagonal. 
Hence, diagonal differences are already incorporated 
in the second order term and we do not have to con¬ 
sider terms of the form D 2 (^i-i,j-i 5 ^i+iy+i, for 
the reduction of anisotropy effects. The model param- 
(3) eters 0 ^ 1 , 0 ^ 3 , 0 ^ 4 , 

A,/^ 2,7 regulate the influence of the different TV 
terms. 

One main reason for considering this anisotropic 
model is that it is computationally feasible for the non¬ 
linear space of combined cyclic and vector space data 
considered in this paper as we will see later on. To 
our knowledge there are no previous algorithms deal¬ 
ing with any kind of second order TV-like problems 
in this nonlinear situation —neither in the isotropic 
nor in the anisotropic setup. As explained, we employ 
diagonal terms to milden unwanted effects caused by 
the anisotropic formulation. 

Next, we consider suitable modifications of the 
above functional to obtain models for the inpainting 
problem with noisy and noiseless data. We start by 
first formulating the inpainting problem. 

Inpainting problem in the presence of noise. 

Given an image domain Uq = {!,..., N} x ,M}, 

an inpainting region U C Uq is a subset of the image 
domain Uq, where the pixel values (i, j) G U, are 
lost. The noiseless or noisy inpainting problems now 
consist of finding a function x defined on Uq from data 
/ given on the complement of the inpainting re¬ 
gion, such that X is a suitable extension to / onto U 
and for the second case additionally denoised. 
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Inpainting without presence of noise. To deal 
with the noiseless situation, we consider the following 
modification of the functional J given by ( 1 ). Since the 
data is assumed to be noiseless, we add the constraint 
that the target variable agrees with the data on the 
complement of the inpainting region. Furthermore, 
the data term considers only those indices for which 
actually data are available. This eliminates the data 
term from the functional. More precisely, the second 
order variational inpainting problem considered in this 
paper reads for a vector space as 

arg min a TV i{x) ^ TV 2 {x) +7 TV 1 1 (x), 

" ’ . (7) 

subject to Xi^j = fij for all (^j) G Or . 

The TV terms TVi, TV 2 , TV 17 defined by (3), (5) 
and ( 6 ) using the difference terms Di, 1 ^ 2 , ^ 1,1 based 
on a norm in the vector space. Due to the constraint 
they actually only act on those difference terms that 
affect an entry in the inpainting region. 

Second order TV formulation of the inpainting 
problem for noisy data. For the inpainting prob¬ 
lem in presence of noise the requirement of equality 
on is replaced by x being a suitable, i.e., smooth 
approximation. In this case, we search for a minimizer 
of the following second order TV functional for cyclic 
data for inpainting: 


for the case of vectors of cyclic and combined data. 
We further unify the notation to D(-; re) for different 
weights w. This overload is employed to avoid addi¬ 
tional notation and should not cause confusion since 
the space under consideration will be clear from the 
context. 

Let w = G w 7 ^ 0, be a vector with 

such a vector w a weight. Spe¬ 
cial cases are the binomial coefficients with alternating 
signs 

For the vectors x^^^ G j = 1,..., d, we employ the 
matrix x := [x^^\ ... to define the abso¬ 

lute finite difference D for these vectors with respect 
to the weight re G by 

d 

D{x;w) = \\xw\\ = I . 

i=i 

For w = bd, we obtain the forward differences of order 
d for the vectors x^^\ ... i.e. 

d-\- 1 

Da{x) := D{x-M = • 

i=i 

Another useful weight used in this paper is re = := 

(— 1 , 1 , 1 ,- 1 ). We denote the corresponding finite ab¬ 
solute difference by Di^ffx) := D{x] 617 ). 


Jq,{x) = FQ^c{x\f) + aT:Yi{x) 

+ ^TV2(a:)+7TVi.i(a:), 

where for any subset 5 C of the image domain, we 
define 


This means we use a data term that enforces similarity 
to the given data / on while applying a regular¬ 
ization based on (3),(5), and ( 6 ) for the whole image 
domain. Specifying the inpainting area as the empty 
set, we obtain the pure denoising problem ( 1 ). 

2.2 Absolute differences for combined cyclic 
and vector space data 

In order to implement the above variational inpaint¬ 
ing problem (7) and the simultaneous inpainting and 
denoising problem ( 8 ) for combined cyclic and vector 
space data, we have to find suitable difference opera¬ 
tors Di, D 2 , Dip for data consisting of combined cyc¬ 
lic and linear space components. In order to do so, 
we first find suitable definitions for vectors of cyclic 
data. Then, we combine these definitions with those 
for the linear space case in a way suitable to the space 
of interest in this paper. We use the symbols D, al¬ 
ready introduced in (3) for the linear space case, also 


Example 2.1. For three points xi^X 2 ^x^ G and 
w = 62 = the second order absolute dif¬ 

ference is given by D 2 {xi,X 2 ,xs) = \\xi — 2 x 2 
cf (4). This can be interpreted as measuring the dis¬ 
tance from y to the midpoint rux •= + xs) of the 

line segment connecting xi and x^; more precisely, we 
have D 2 {xi,X 2 ^xff} = 2||-h xffj — ^ 2 1|. The situa¬ 
tion is shown in Figure 1 forn = 2, but also illustrates 
the situation for general n > 2 in the plane defined 
by xi,X 2 ^xs. For n = 1 the situation simplifies to X 2 
always lying on the line —though not necessarily the 
segment — connecting xi and X 3 . 

We consider the cyclic case next. Let denote the 
unit circle := {pi -\-P 2 = I : p = (pi,P 2 )^ G 
endowed with the geodesic or arc length distance 

(P, q) = arccos (p, g), p, g G 


Given a base point g G , the exponential map exp^: 
M ^ from the tangent space TqE>^ M of at g 
onto is defined by 


expg(x) = RxQ, 


( cosx 
ysinx 


— sinx 
cosx 


This map is 27r-periodic, i.e., expg(x) = expg((x) 27 r) 
for any x G M, where (x) 27 r is the unique point such 
that 


X = 27 r/c + (x) 27 r with (x) 27 r G [ —tt, tt), /c G Z. (10) 
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Figure 1. The three points xi,X2,X3 G illustrate the 
multivariate finite difference D2{xi, X2, xs) = 2||mcc — 
X2II, i.e. the distance of X2 to the midpoint rux 
|(xi + X3) of xi and X3. This measures how “near” 
they are to lying equally distributed on a line segment 
in the right order. 

If we fix g, we obtain a representation system of 
i.e., expg is a bijective map where expg(O) = q and 
there is a unique x G [—tt, tt) for each p G such 
that expg(x) = p. A vector x G [—7r,7r)’^ represents a 
point in by component-wise application, and for 

a point q G the map exp^: [—7r,7r)’^ ^ 

expq{a;) := (exp^^(a;i),... ,expq^(a;„j))'^, is bijective 
where the properties from above hold component-wise. 
A distance measure on is given by 

d(si)^{p,q) = ||(arccos(pi,gi))™ J|. 

In the following, we introduce higher order differ¬ 
ences on We employ the representation system 

= [— TT, tt) induced by using an arbitrary but fixed 
exponential map. Let x^^^ G [—7r,7r)’^, j = l,...,(i. 
Using the notation x := {x^^\ ..., x^^^) G [—tt, 
the absolute eyelie differenee of x^^\ ... with re¬ 
spect to a weight w G ]R^\{0} is defined as 

D{x] w) := min D{[{x^^^ -h (a,..., x^^^ -h <a)] 27 r; 

aER'^ 

where [y] 27 T foi* some ^ G is multivalued and its 

ith component ([^] 27 r)i is given by 

(^i) 27 r, if Vi 7 ^ + l)7r for all z G Z, 

Ttt else. 

( 11 ) 

This definition may seem a bit technical at first glance. 
However, it allows for two points x^^\ x^^\ i^j G := 
{1,... ,(i}, having the same value x[^^ = in one 
component I G {!,..., m} to be treated differently, cf. 
the definition for in [6, Section 2]. In fact, we may 
choose any g G as a base point for our represen¬ 

tation system, which shifts any set of points given with 
respect to exp^/ by a fixed value of a := exp“^(g'). 
When the shift by a G is small enough, such that 
no component of x is affected by the component-wise 
application of [•] 27 r, both representation systems yield 
the same value. Using this notation we can simplify 
both the definition of the second order difference and 


the proximal mappings derived later on. The mini¬ 
mum in the definition of the difference simplifies to 

D{x]w) = min D{[{x^^^-Xk-^ir ,... ,x^^^-Xk^7r)]27r]w), 
kei^ 

where Xk •= This is illustrated in Figure 2 

for three points x^y^z G (S^)^. We note that this 
definition contains the notion introduced in [6] for 
as a special case. 

Finally, we come to the space of interest in this pa¬ 
per which is X := x In this space a vec¬ 

tor X = {xi)^^^ consists of two parts: the phase¬ 
valued xg := {pi)^i G and the real valued 

components xr = ^ of x e X. In 

the following, we will use any representation system 
in order to write the cyclic components xg as a vector 
xg G [—TT, ir)^. This also allows for x G A to be seen as 
a vector, where the first m components are restricted 
to [— 7 r, 7 r) and the remaining n ones are real-valued. 
The distance of two points x,y G X on this product 
space is given by 

dx{x,y) = ^J\\xw-m\\'^ + d(si)^ixs,ysr. 

For a set of points x^^\ ... G A, using the nota¬ 
tion X = {x^^\ ..., x^^^) as before, the finite differenee 
for eyelie and noneyelie data with respect to a weight 
w G ]R^\{0} is defined by 


D{x;w) := ^/D{xJ^;w)^ + D{xg;w)^. 

We further introduce the short hand notations 

Dd{x) = D{x, hd), X G d G N, (12) 

to denote the corresponding absolute finite differences 
of order d. Furthermore we introduce —with a slight 
abuse of the difference notation— the seeond order 
mixed differenee of four points, e.g. given on a 2 x 2 
subset of the pixel grid ( 8 ) 0 ) by Di^i{x) := D{x;bi^i) for 
X G with = (—1,1,1,—!)^. For the weights 
corresponding to first and second order differences, we 
have a particularly nice representation, which is given 
by the following Lemma. 

Lemma 2.2. For w G {^1,^2, ^1,1} ci'^d x G X^ 
where d denotes the length of w, we have 

D{x,wf = ||(a;sto)27r|P + (13) 

Proof. For the real-valued components there is noth¬ 
ing to show. Hence we may restrict to n = 0 (which 
corresponds to purely cyclic data) and thus have to 
show that D{x^w) = ||(xre) 27 r|| for x G (S^)’^^^. To 
this end we apply Proposition 2.5 of [ 6 ] to each row x 
and conclude the validity of (13). □ 
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(a) Different shifts, where the first two, Fs,Fs' yield the (b) Three representation systems with points x,y,z at 
same constellation of x^y,z. ±7r in at least one dimension. 

Figure 2. For given points x,y,z G [—tt, tt)^ in a representation system Fs, i.e., with base point s the two subfigures 
illustrate different shifts: ((a)) shifts by arbitrary a G e.g., to s' = s + a yielding the same value of D 2 , 
and ((b)) cases yielding different values for the second order difference and a minimum occurring for x,y, z at the 
borders of the representation system. 


2.3 Inpainting combined cyclic and vector 
space data 

We can now apply the definition of absolute differ¬ 
ences for combined cyclic and vector space data we 
derived in Section 2.2 to obtain variational models for 
inpainting and simultaneous inpainting and denoising. 
This extends the models for the Euclidean differences 
in Subsection 2.1 using a first order TV term (3) as 
well as to the second order TV terms (5) and (6) to a 
more general setting. The noiseless inpainting model 
now reads 

argmin (aTVi(x) + f3Ty2{x) + 7TVi,i(x), 

subject to Xi^j = fij for all (i, j) G 11^. 

(14) 

Here, TVi is defined by (3) incorporating the first or¬ 
der absolute differences Di given by (12) and the sec¬ 
ond order TV terms TV2, TV 1,1 are defined by (5) 
and (6) employing the second order absolute differ¬ 
ences D 2 , Dl l given by (12), respectively. 

Proceeding similarly, we get a variational formula¬ 
tion of the inpainting problem for noisy combined cyc¬ 
lic and vector space data by computing a minimizer 
of 

Jn{x) = FQc{x;f) + aTYi{x) 

+ /3TV2(:r)+7TVi,i(a;). ^ ^ 

Here the data term is given by (9) using the distance 
function on V = (S^)"^ x As in the noiseless sit¬ 
uation, TVi is defined by (3) again incorporating the 
first order absolute cyclic differences Di and the sec¬ 
ond order TV terms TV2, TVi^i are defined by (5) 


and (6) employing the second order absolute cyclic 
differences D 2 , Di^i from (12), respectively. 

3 Algorithms for inpainting and 

denoising combined cyclic and linear 
space data 

In this section, we derive algorithms to solve the in¬ 
painting problem (14), and the combined inpainting 
and denoising problem (15). Note that the latter in¬ 
cludes the denoising of combined cyclic and vector 
space data for the case of H = 0, an empty inpainting 
set. These algorithms are based on a cyclic proximal 
point algorithm whose concept we recall in Section 3.1. 
We derive explicit formulas for the proximal mappings 
that are needed for inpainting and denoising of such 
combined data in Section 3.2. Using these explicit 
representations, we derive a cyclic proximal point al¬ 
gorithm for inpainting noiseless combined cyclic and 
vector space data and similarly for simultaneously in¬ 
painting and denoising data in Section 3.3. This also 
includes an efficient choice for the cycles involved. Fi¬ 
nally in Section 3.4, we prove convergence of our algo¬ 
rithm to a minimizer under certain conditions that re¬ 
flect the space-inherent non-convexity of the involved 
functionals. 

3.1 The cyclic proximal point algorithm 

For a closed, convex and proper functional ip: ^ 

M U {00} the proximal mapping is given by 

Pi'OX;,^(/) := argmin 1||/ - x\f‘ + \^{x), f e R”, 
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where A > 0 is a tradeoff or regularization parame¬ 
ter. The fixed points of pTox^^{f) are minimizers of 
(fi. Hence, if the proximal mapping prox;^^(/) can be 
computed in closed form, an algorithm for finding a 
minimizer is given by iterating 

for some starting value This algorithm is 

called proximal point algorithm (PPA) and was intro¬ 
duced by Rockafellar [68]. It was recently extended 
to Riemannian manifolds [35] and also to Hadamard 
spaces [3]. Denoting the distance on a Riemannian 
manifold M. by (Im the proximal mapping on a man¬ 
ifold reads for a function (f\ A4^ ^ M as 

Prox^^(f) := argmintd_M(f,xf + A(p(x), f € M"". 

xeM^ z 

We note that on some manifolds, e.g. the spheres 
there is no definition of a (globally) convex function. 
Hence the minimizer might not be unique. This is for 
example the case when looking at our space as 
this is not even the case for cyclic data, cf. [6]. 

If the function (p can be split into simpler parts, i.e. 
p = for which then individually the proxi¬ 

mal mappings are known in closed form, a similar al¬ 
gorithm is given for a sequence {Xk}k of regularization 
parameters by 

l = k = l,2,..., 

and it is called cyclic proximal point algorithm 
(CPPA). Its formulation on Euclidean space is derived 
in [9], see also the survey [8]. It converges to a mini¬ 
mizer of p if 

E OO r —^ OO 

k^o^k = oo, and 2^^^Pk<oo. (16) 

The concept of CPPAs for Hadamard spaces has been 
treated in in [2]. A CPPA for TV minimization for 
manifolds and in Hadamard spaces has been derived 
in [85]. For second order TV type problems, a CPPA 
to denoise data was derived in [6]. A preliminary 
model, different to the one appearing in this paper, 
was applied to inpainting of data in [7]. For man¬ 
ifold data in general, the main challenge is to derive 
proximal mappings which are as explicit as possible. 

3.2 Proximal mappings for inpainting 

Here we derive closed form expressions for the prox¬ 
imal mappings needed to make the cyclic proximal 
point algorithm from Section 3.1 work for the inpaint¬ 
ing problems (14) and (15) in the nonlinear spaces con¬ 
sidered in this paper. In particular, we derive proximal 
mappings incorporating constraints directly. 


We first need some basic results on the linear case 
which involves vectors of real-valued data only. To 
this end, we start with a generalization of [6, Lemma 
3.1]. We derive explicit expressions for the proximal 
mappings of functions living on linear spaces which 
are of the form 


where the target variable x is a matrix in and 

where d corresponds to the length of w. The vector a 
introduces an offset. We employ the notation ||^||f = 

denote the Frobenius norm of a 

matrix y. 

Lemma 3.1. Let f = (/^^\ ...,/^^^) G be a 

matrix whose columns represent the data vectors, 
lett) ^ w (w not necessarily a weight), and A > 0 
he given. For the functional 


E{x]f,a,w 


|/-a;|||. + A||a;w-a||, (17) 


with target variable x G the minimizer x is given 

by 

X = f — msw^, (18) 


ll/-«ll ^ 0, 

0 else, 


and m := min{A, 
E{x; /, a, w) is given by 

E{x]f,a,w) 


Ih||2(p2 + A(| 


The minimum 

(19) 

if m < X, 

— A)) otherwise. 


Furthermore, given data /, / G and different 

offsets a,d e the following implication holds: 

\\fw-a\\<\\fw-d\\ (20) 

min E{x;f,a,w)< min E{x;f,d,w). 


Proof. We first reduce the functional to be minimized 
to an equivalent problem without offset. By assump¬ 
tion there is an index j such that Wj 7^ 0, which allows 
us to write (17) as 

E{x-,f,a,w) = h\f - x\\l + X\wj\\\{x - f-ej){^)\\ 

UJ j UJ j 

Defining the new target matrix y := x — :^ej, the 
new data matrix g = / —the new regularizing 
parameter v := AlrCjH, and the new (not necessarily 
weight) vector v we obtain the new problem 

F{y,g,v) = t\\g-y\\l + p\\yv\\, (21) 


8 






where the second term is free of an offset. The relation 
between minimizers x of E and ^ of F is given via 
y = X — —ef. 

We now consider the problem ( 21 ) and first 
show (18) for F. The corresponding statement for E 
follows by carrying out the resubstitution. If = 0 
then we have F{g; g^v) =0 and hence ^ ^ is the 

minimizer of ( 21 ). So we may assume liar’ll 7 ^ 0 in the 
following. We now distinguish whether liar’ll 7 ^ 0 or 
liar’ll =0. In the first case, we may differentiate F, 
and setting the gradient of F to zero results in 


f^ = y-9 


V 

w 


{yv)v^. 


We multiply by v to obtain yv — gv = —||i;|p. 

Rearranging yields (l = gv^ which implies 

i.e., both vectors have the same direc- 

\\yv\\ \\ 9 v\\ ’ ’ 

tion. 

This leads to 


y = g-iy 


gv 

Il5v|| 


V 


T 


For liar’ll = 0 we look at the subgradient of F. As 
condition for a minimizer y, we have that ^ ^ is in 

the subgradient of For y with \\yv\\ = 0, this 

subgradient is given as {zv^: ||z|| < z/}. When con¬ 
sidering the functional F, the amplitude m from the 
assertion of the lemma reads m = min{z/, ||5''^||/||'^|| }• 
If 115 HI/ 11 HI < z/, then F is differentiable at ^ = 
g — msv^, 11 ^ 1 ; 11 7 ^ 0 , and we are in the previously 
considered case. Hence, we may assume that m = u. 
Then, for ^ ^ — msv^, we have y — g = msv^ with 

m = V. This shows that y fulfills the condition of a 
minimizer. In consequence, (18) is true for the func¬ 
tional F. Then resubstituting shows (18) for F, and 
plugging X into E we get (19). 

It remains to show the implication (20). To this 
end, let g = and ji := By the assumption 

of (20), \\g\\ < ||/2||. We consider three cases. If \\jl\\ < 
A, then \\g\\ < A. Hence, the minimizer of E{x; /, a, w) 
equals ^ H'a^lli IImIP? and the one of E{x; /, d, w) equals 
iMMyr > hMUyf which shows (20). If m > 
A and ||/i|| < A, we have to consider the second line 
of (19) for ||/i||. From this we obtain a minimal value 
of E{x; /, d, w). We have to show that 




but this is a consequence of the second summand on 
the left hand side being positive and A > \\g\\. Finally, 
if both \\fl\\ > A and ||/r|| > A, we apply the second 
line of (19) and see that we need ||yd|| — A > ||/i|| — A 
for the statement to hold. This is true by assumption 
which completes the proof. □ 


Example 3.2. We continue with the situation from 
Example 2.1 and take three points G j = 1, 2, 3 


and denote f = (/^^\/^^\/^^^). Depending on the 
chosen value for A in the proximal mapping, there are 
two possibilities: If m = < A, we obtain 

three points x = (/) = / — msb^ that lie on a 

line, cf Figure 3 (a). If m > X, then the result x of 
the proximal mapping does not yield D 2 {x) = 0, but 
the 'movement of the points in direction s is restricted 
by A 6 j, cf Figure 3 (b). 

After these preparations, we now deal with the prox¬ 
imal mappings needed for the inpainting problems (14) 
and (15) for combined cyclic and vector space data. 
In particular, each data item now is an element of 
A = (S^)"^ X As motivation, let us first have a 
look at the first order difference Di and the inpaint¬ 
ing problem (14) for noiseless data. By the constraint 
Xij = fij outside the inpainting region, it might hap¬ 
pen that at the boundary of the inpainting region the 
member Xij = fij is fixed but its neigbor, say Xi^j 
may vary. Then, we have to study the correspond¬ 
ing functional Di{xij,Xi^jj^i) for fixed Xi^j and find its 
proximal mapping. The following theorem deals with 
this issue in a more general setup. 

Before we state the theorem we introduce some 
notation: By again using a representation system 
= [—7r,7r), we can interpret any x G A as a vec¬ 
tor X = {xg,XM.)^ G where the same representa¬ 
tion system is used for all cyclic components. 

We define (•)a’ : ]R(^+^)x^ ^ by 

{x)x := ((^s)27r,^R)^, 

where (xs) 27 r is defined as the component-wise appli¬ 
cation of (•) 27 r as given in ( 10 ), i.e. is applied to the 
m cyclic components of each column vector x^^^ G A. 
Similarly we define 


[x]x := ([a;s]27r,a;R)'’’, 

where [•] 27 r defined in ( 11 ) is applied to the phase¬ 
valued components of x. 

In the following we consider a weight vector re G 
a data matrix x = {x^^\ ... with each member 

x^^^ having values in A = x and a subset 

A C {1,... ,d}. We partition w into a variable part 
and into a fixed part w according to whether the 
index i of Wi belongs to A or not. Accordingly, we 
partition x into a variable part x^ and into a fixed 
part X and consider the mappings 

ipA‘ ^ D{x, w) ( 22 ) 

for the corresponding differences D{x,w), where only 
the x^ are considered as variable and the x are fixed 
values. We derive an explicit representation for the 
corresponding proximal mappings in the following the¬ 
orem. 
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(a) Minimizing D2{f) = 2 ||m/ — to (b) Tradeoff between minimizing D2 and a maximal 

zero by moving / to x, where = rux. movement of /b\ j = 1, 2 , 3 . 


Figure 3. Two cases of the second order difference D2{f) in by looking at the plane generated by G 

R^: The original points / are moved onto x towards forming a line: ((a)) yielding D2(x) = 0 , i.e. X2 is the mid 
point rux of xi and X3, and ((b)) reducing D2{x) < D2{f) but restricting the movement to be less then A G R 
for and 2 A for respectively. 


Theorem 3.3. Let w be one of the weights w = 
(—1,1), w = (1,—2,1), or w = (—1,1,1 ,—!) whieh 
eorresponds to eonsidering the first order differenee 
Di and the seeond order differenees D 2 and Di^i, re- 
speetively. Let d be the respeetive length of w, A C 
{ 1 ,... ,(i}, and w be partitioned into the eorrespond- 
ing variable part and into a fixed part w. We par¬ 
tition X, / G aeeordingly and let f = x. Then, the 
proximal mapping of if a defined in ( 22 ) is given by 

(/“) = {r -ms {w°-A)x, 

with the parameter A > 0 ; here, the direetion(s) s and 
amplitude m are given by 

wlfZlxW 

0 else, 

and 

m = mini (23) 

Remark. We note that the braeket [•];f and thus the 
proximal mapping (having an additional value for eaeh 
additional instanee of s), is multivalued if some eom- 
ponents of {fsw) 27 r are equal to —tt. More preeisely, 
if there are I G {1,..., n} sueh eomponents, we obtain 
if solutions from the different instanees the veetor s 
might take. The reason for this is, that the mapping 
if A is no longer eonvex for data in X, and that the 
minimizer defining the proximal mapping is not neees- 
sarily unique. Owing to this observation, we eonsider 
set valued proximal mappings gathering all minimiz- 
ers. We notiee that the above proximal mapping is 
single-valued if and only if {fgw) 27 r ^ (— 7 r, 7 r)^. This 
is the generie ease. The degenerate ease involving an¬ 
tipodal points appears rather seldom in praetiee; at 
least, in a non-artifieial noisy setup, it is very un¬ 
likely to eneounter antipodal points. Then, at least the 


proximal mapping is single-valued yielding a determin- 
istie result. However, we notiee that jumps of hight 
elose to TT are problematie sinee then stability issues 
appear. Furthermore, data with antipodal points or al¬ 
most antipodal points may often be interpreted as not 
fine enough sampled data. This means, if the sampling 
rate is higher, the distanee of nearby data items might 
get smaller and the situation might just disappear. We 
note that this does not exelude the possibility of jumps 
in the finer sampled data. The point is that large erit- 
ieal jumps might be revealed as smaller jumps. 

Proof of Theorem 3.3. In order to derive explicit for¬ 
mulae for the proximal mappings, we have to find the 
minimizer(s) of 

r ,«,):= 1 ^ + XDix-, w). 

jeA 

By Lemma 2.2 we may rewrite Ex{x^-, f^yw) as 

£xix^-,r,w) 



+ min Ay||a;Rw ;||2 + \\xsw - 2Trap. 

We now use that x = f and employ the notation H-Hf 
for the Frobenius norm. For the remaining part of the 
proof, let t := \A\. We obtain that 

= min (xII/k-^IIf + xII/s -2;s-27rfc||F 

+ + /rw||2 + - {2TTa - /sw)||2 
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Using this, we rewrite the minimization problem to 
min ,w) 

a:“e[-7r,7r)’^x*xR’^x* 

= min min Ek(j{x^;f,w) 

aEZ'^ 


= min 

kez'^^'^ 

(jEZ'^ 


min 

cc«e[-7r,7r]^xtxR^xt 


Ek,a{x°'-,f,w)l24) 


i G {!,..., m}. Then ||z^/c,a|| is minimal with respect to 
k, a if and only if kw^Ecr = r where r = (ri,..., r^)^, 
cf. (27). For each fixed matrix k G there is a 

uniquely determined vector a = d-{k) solving this sys¬ 
tem of linear equations. Such a pair minimizes ||z^/c,cr|| 
w.r.t. /c, a and 


= (Mx- 


where 

Ek,a{x"‘-, f,w) 

= 2 11-^* ~ ^rIIf + 2^\fi ~ ~ 27r/c|||- 

+ + /k^P + \\x^w<^ - (27rcr - /sw)P. 

(25) 

Having a look at (25), there exist values k, a such that 
a(^) ~ ^k a(^)' Summing up, the problem reduces 
to finding the minimizers of all Ek^a in [—tt, and 

comparing their value. For the remaining part of the 
proof, let On be a zero (column) vector of length n and 
On,t be a zero-matrix of dimension n x t. 

For any k, a, the functional E^^a- has a unique min- 
imizer given by Lemma 3.1 as 

xl,c- = /“ - 27r (yQ ^ - Sk,amk,a{w°')'^. (26) 

We first derive Sk,a using Lemma 3.1, where we use 
the notation, 

^k,a — ^/c,cr/IIII 2 • 

We notice that the data for Lemma 3.1 is given by 
— 27r ( 0 ^ ^) and the offset a in the same lemma is 
a = 27r ( 0 ^ ) - fw. We get 

n,a = (/“ - 27 r (0^ )) + /tZ> - 27 r () 

= fw- 2 TT (( 0^ ) + ( o"„ )) , 

and 

{ Uk 

TTZ; 

By (20), we have to find the minimum of the norms 
||z//c,cr|| with respect to k^a in order to find the mini¬ 
mum or minima and their corresponding minimizers of 
the Ek^cj‘ We first consider the following special case, 
where 

{fw)i ^ 27rZ — TT for all i G {1,..., m}. (28) 


Using (26) we obtain the corresponding minimizer 
w.r.t. X as 

*l,, = r- 2 » 

with m, 5 as given in (23). Now there is precisely one 
/c* with its corresponding a = (j(/c*) such that x^^ ^ G 
[— 7 r, 7 r)’^^^ X which implies that 

X* =xt,^^ = {r-ms{w^f)x 

is a minimizer of £x by (24). Concerning uniqueness 
we notice that a i-G ||z^/c,cr|| has a unique minimizer. 
By (20) this implies that one may minimize w.r.t. /c, x 
in (24) choosing a = a{k) as previously in this proof. 
This unique minimizer is x* since it is the only can¬ 
didate with its first m components in [— 7 r, 7 r). This 
finishes the special case (28). 

For the general case, let G C {1,..., m} be the set 
such that 

{fw)i G 27rZ — TT for i G G 

If G is empty, we are in the case (28) we considered. 
Furthermore let G^ := {1 ,..., m}\G. Hence for each 
i G G^ we find Vi such that {fw)i — 27rri G (— 7 r, 7 r) 
with the same arguments as for the case (28). For 
each i G G there is G Z such that {fw)i = (2r^ — 
l) 7 r. Then attains its smallest value exactly when 
kiW^ E (Ji = Vi, for all i G and when kiW^ + G 
— 1}, foi* all i e G. Following the same steps 
as above, let us fix k and determine all au = d-{k,u), 
u G {0, l}^,Ui = 0 for i G G^, that fulfill one of the 
systems of equations 

cr(^, u) = r — kw^ — 14, where r = (ri,..., r^)^. 

These are 21^1 systems systems of equations, one 
for each value of i 4 , each having a unique solution. 
By (20), we only have to consider the functionals E^^cr 
corresponding to such a pair of parameters (k^dy)- 
We get that the components of ^ with indices in 
G^ are given as in (27), whereas the components with 
indices in G are given by 



More precisely, in none of the cyclic data dimensions 
1 = 1,..., m of / G the scalar product with the 
weight {fw)i is an odd multiple of tt. There exist 
ri, ..., G Z such that {fw)i — 27 rr^ G (— tt, tt) for all 




i+U 


From ^ , we obtain 5 ^ ^ again by using (23) for this 
special z/. We furthermore notice that ^ equals 
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the m given in (23) and is especially independent of 
the particular choice of u. Hence we get 

= /“ - 27r (J^ 1 - St - 

Now we follow the lines of case (28) to conclude that 
where 5 is one of the 2 l^l in¬ 
stances of ; and these are precisely the mini- 

mizers of £{x^] ^w) w.r.t. lying in [— 7 r, 7 r)"^^^ x 

I^nxt^ This completes the proof. □ 

Example 3.4. ITe continue with the situation from 
Examples 2.1 and 3.2, but split the three points f = 
into an active and fixed part using 
A = {2,3}. In other words f^^^ is fixed, we have 
^(1) = j(i)^ and both fA3) and fA,^) are active, i.e. 
^( 2 ) ^^(3) affected by the restricted proximal map. 
The movement again depends on the chosen value for 
A which is also illustrated in Figure 

We also need proximal mappings for the data term 

Fix;f). 

Proposition 3.5. For f,g^ k = NM, let 

k 

i=l 


A splitting for noiseless inpainting. Each sum¬ 
mand in the first and second order differences in (14) 
is incorporated into a proximal mapping using Theo¬ 
rem 3.3 by setting the values affected by the subjection 
as fixed and keeping the remaining ones as active. 

If two summands act on distinct data, their proxi¬ 
mal mappings can be computed in parallel. This re¬ 
duces the cycle length c of the CPPA tremendously 
and provides an efficient, parallel implementation. In 
the following, we will split each of the summands 

J{x) = oTVi(x) +^TV 2 (x) + 7 TVi,i(x) 

into 


where the summands for the first order horizontal dif¬ 
ference are split as 

Jl + J 2 •= OLi ^ Di{x2iJ,X2i^lj) 

ihj) 

T l(^ 2 i+l,ji‘ 5 ^2i+2,ji‘)5 

(hi) 

which similarly yields J 3 , J 4 for the vertical first order 
difference and J 5 ,..., Jg for the two diagonal sums. 
The second order term is decomposed into three sums 
J 9 , Jio, Jii, given by 


Then, the minimizer(s) x of £ are given by 
^gTXf ^ A 




+ A 1 + A 


2ttv 


),■ 


where v = 


k,n-\-m 

hi=i 


G is defined by 


I 


0 

0 


sgn{gP 



if j > m, 

if 32 -/bl < m, 

if\g2 ~ 4*^1 > < m. 


J 9 '■= /3l y] D2{x3i-ij,X3ij,X3ixij), 

(hi) 

-^10 •= A ^ ^ 2 (^ 82 ,i, ^3i+l,i, ^3i+2y), 
(hi) 

Jll (^1 ^ ^ D 2 (^3i-|-l,i ; ^3i-|-2,jf; ^3i-|-3,i ) ; 

(hi) 


and analogously we obtain J 12 ,..., J 14 for the vertical 
second order difference. Splitting the term 7 TVip(x) 
in this manner yields J 15 ,..., Jig as follows: 


Proof. We observe that the functional £ under consid¬ 
eration only involves squared distances. This can be 
handled separately for every index {i,j) G {1,... ,k} x 
{I,..., n + m}. Hence, the proposition follows from 
considering both the M-valued data and the -valued 
data case separately. This has been done in [ 6 ] in 
Lemma 3.3 and Proposition 3.7, respectively. □ 

3.3 Cyclic proximal point algorithms 

The proximal mappings from Theorem 3.3 can be 
efficiently applied in parallel to compute minimizers 
of (14) and (15) using a cycle length in the proximal 
point algorithm from Section 3.1 by splitting the func¬ 
tionals accordingly. 


Ji5 = 7 E^i, l(^2i,2i, ^2i+l,2j, ^2i,2j+l, ^2i+l,2i+l), 

(hi) 

J16 = 7 E^. l(^2i+l,2j5 ^2i+2,2j5 ^2i+l,2j+l5 ^2i+2,2i+l)5 

(hi) 

Jl7 = 1 ^l,l(^2i,2j+l5 ^2i+l,2j+l5 ^2i,2j+2, ^2i+l,2j+2)5 

(hi) 

Jis, = 7 -Pl,l(^2i+l,2j+l, ^2i+2,2j+l, 

(hi) . 

^ 2 i+l, 2 j +25 ^ 2 i+ 2 , 2 i+ 2 )• 

In any of these functions J/, / G {I,..., 18}, any data 
point Xi^j G T is occurring at most once and hence all 
proximal maps a function Ji consists of can be evalu¬ 
ated in parallel. This leads to a cycle length of c = 18. 
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j(a,2) j(a,2) 




(a) Minimizing \D2 < \ with fixed (b) Minimizing |D2 > A with fixed 

Figure 4. Minimizing a second order difference in where is fixed, i.e. and A = {2, 3} are the active 

data points. Again, two different values of A are shown: ((a)) A > |D2(/;ie), i.e. the corresponding proximal 
mapping reaches the minimum, ((b)) A < |D2(/;re), i.e. the corresponding proximal mapping is just a step 
reducing the value D2{x;w) < D2{f;w). 


A splitting for combined inpainting and denois- 

ing. In order to derive a cyclic proximal point al¬ 
gorithm for the combined inpainting and denoising 
model (14), we encounter two differences compared 
to the previous derivation: All data is always marked 
as active because no index (i,jf) G is restricted by 
a constraint and we further obtain a data term, i.e. 
we additionally have 

( 29 ) 

which can be evaluated in parallel using the proximal 
mapping given by Proposition 3.5. 

Initialization. In order to initialize the algorithm, 
we employ an the idea of unknown boundary used 
in [1], which can easily be implemented during the 
first iterations of the CPPA: all unknown values 
Xij, {i^j) G ft are initialized whenever setting a corre¬ 
sponding difference D{x,b) = 0 yields a unique value 
for this pixel. Afterwards, this data item is set to be 
known and can be used to initialize other unknown 
pixel. Hence after at most k = max{ A, M} iterations, 
all pixels are known. 

The complete procedure for both models of noiseless 
and noisy inpainting is summarized in Algorithm 1. 

3.4 Convergence Analysis 

As typical when dealing with nonlinear geometries, 
we show convergence under certain conditions. These 
conditions are comparable with the ones employed 
in [6]. In particular, the data is assumed to be dense 
enough in the sense quantified later on. This means 
that the data has to be locally nearby which does not 
mean that circular data components are (globally) re¬ 
stricted to certain sectors - the data in these com¬ 
ponents my wrap around. Similar restrictions on the 
nearness of data and even more severe restrictions re¬ 
quiring almost equidistant-data have been imposed in 
the analysis of nonlinear subdivision schemes; see, e.g.. 


[46, 82, 84]. As it is also pointed out in these references, 
the analysis is qualitative in the sense that empirically 
convergence is observed for a significantly wider range 
of input data. 

Compared with the pure denoising setup considered 
in [6], there are several issues we have to deal with in 
the inpainting situation here: first, the proof there 
relies on the uniqueness of the minimizers in the un¬ 
wrapped situation which is not given for inpainting; 
second, a main step in the aforementioned proof is 
based on bounding the distances dx{fi^j^Xi^j) for all 
pixels to get information on x, while for inpaint¬ 
ing the values fij^ (bi) ^ missing. In addition 

we consider a more general data space. 

We first discuss the conditions we impose for our 
convergence analysis. Then we derive the necessary 
information to prove our main result formulated as 
Theorem 3.14. It states that both the algorithm pro¬ 
posed for inpainting and the algorithm proposed for 
simultaneous inpainting and denoising converge to a 
minimizer. 

We employ the following notation to denote the dis¬ 
tance on the first m components of two data items 
x^y G A. We notice that those are the -valued com¬ 
ponents. We let 

dx,m{x,y) ■■= d(si)m{xs,ys). 

Our first condition is that the data /; O'" —>■ X given 
on the complement of the inpainting region ft is dense 
enough in the sense that the distance between pixels 
and their neighbors in ft^ is sufficiently small (in the 
respective spherical components). 

In order to give a precise definition of this we re¬ 
quire some preparation first. On the domain grid, 
we consider the distance d{{i, j), (kj)) which is the 
length of the shortest path with respect to the eight- 
neighborhood (with diagonal distance y/2) connecting 
the indices (i, j) and (/c,/). We consider a covering of 
the image domain r^o with balls {H((i, j), r^^j) : (i,j) G 
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Algorithm 1 CPPA for minimizing (1) for x valued data 

Input a sequence {Xk}k of positive values, cf. (16), 
parameters a = ( 01 , 02 , 03 , 04 ), 13 = {Pi, P 2 ), 1, 
a set C Ho? and data / G 

function CPPA(o, p, 7, {Xk}k, f) 

Initialize = fij, {Pj) ^ as active, as unknown for {i,j) G H and k = 0 

Initialize the cycle length as c = 18 (noiseless case) or c = 19 (noisy case) 

repeat 

for / 4 — 1 to c do 

^ prOX;)^^ (x^^+^^) 

employing Theorem 3.3 and Proposition 3.5 

k ^ k 1 

until a convergence criterion are reached 

return 


centered at pixels (i,jf) G in the complement 
of the inpainting area. Ne note that the radius of 
the ball B{{i, may vary with {i,j) G and 

that balls are formed w.r.t. the distance d{{i,j), {k, /)) 
introduced right before. We require that the graph 
induced by this covering is connected. Here the in¬ 
duced graph is formed by using the indices in as 
vertices and by connecting to vertices {i,j) and {k,l) 
whenever {k,l) G B{{i,j),rij) or vice versa meaning 
that (i, j) G B{{k,l),rk,i). 

Using such a covering {H((i, j),r^j) : (i,j) G 
as well as the shorter notation Nij = B{{i, j),rij), we 
let 

' (.SS- A-'*- 

We emphasize that imposing a bound on (30) still al¬ 
lows for jumps which are not too high. As already 
pointed out, such restrictions on the nearness of data 
are typical for the analysis of algorithms of nonlin¬ 
ear data in general; cf. [6, 46, 82, 84]. We note that 
d^{f) depends on the chosen covering and that we 
suppress this dependence in the notation. The above 
definition takes the inpainting region into account and 
only restricts the spherical components. It turns out 
that for the non-spherical linear space components, no 
restrictions are necessary, and that large distances in 
these components do not influence the behavior in the 
spherical components negatively. When there is no in¬ 
painting region, i.e. U = 0, and we are in the pure 
denoising situation, we fix the covering by fixing Ni^j 
to be the eight-neighbourhood of (i,j) and use the 
notation 

doo{f) = d^if). 

Using this notation, we assume that (i) the quantity 
d^{f) is small enough; precise bounds on d^{f) are 
given in the lemmas and theorems later on. 

Our second requirement is that (ii) the parame¬ 
ters a,p,j are sufficiently small. For large parame¬ 
ters, solutions become almost constant which is often 


undesired and causes an interpretation problem, e.g. 
when the original data is equally distributed around 
the circle. Finally, we require (iii) that the parameter 
sequence {Xk}k of the CPPA fulfills (16) with a small 
norm. The latter can be achieved by rescaling the 
parameter sequence. 

Our analysis is based on an unwrapping procedure 
which means that we ‘lift the whole setup’ to the uni¬ 
versal covering of A which we denote hy y = 

Similar to the notation we use 

k - y\y,m ■■= WixV'jLi - (%)7ill 

to denote the distance on the first m components of 
the data x,y ey. 

Universal coverings stem from algebraic topology. 
We refer to [47] for an introduction. A covering con¬ 
sists of a covering space and a canonical projection 
(inducing discrete fibers). We here explicitly consider 
the canonical projections tTx which are for x G A given 
by 

■^x{y) = ■^x{ys,yR) = {exp^s{y§),yR), y ey, 

i.e. the linear space components remain unchanged 
and the cyclic components undergo the exponential 
mapping component-wise meaning that 

exp^s(ys) = (exp^^yi),... ,exp^^(y^)). 

It is well known that continuous mappings to the base 
space have a lifting to the covering space. The lifting 
is uniquely determined by specifying 7r~^{x) for only 
one X. This lifting construction also applies to discrete 
mappings g: Qq ^ At defined on the rectangular grid 
Uq whenever doo{g) < tt. We record this observation 
for further use. 

Lemma 3.6. Let g: Lto At be an image 

with doo{g) < TT, and eonsider g G A fulfill¬ 
ing dgi{qi, {gi^i)i) < IT, i = l,...,m, i.e. no pair 
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of cyclic data components is antipodal. We choose 
gi,i ^ y such that 7rq{gi^i) = gip. Then there exists 
a unique lifted image g\ TLq ^ y such that 7Vqg = g 
holds component-wise and doo{g) < tt. 

Next, we lift the inpainting functionals and derive 
relations between the lifted and not lifted functionals 
and lifted and non lifted discrete functions. To pre¬ 
cisely formulate these relations we need some prepa¬ 
ration. 

For ^ > 0 and data /: X given on the com¬ 

plement of the inpainting region f], we consider the 
class S^{f^ 5) of grid functions x defined on the whole 
domain which we define by 

S^{f,5) = {x-.nQ^X-. eoo{xJ) < 5 }, 


Lemma 3.7. Let f: X with d^{f) < be 

given and let q G X be a point not antipodal to fv{i,i) 
in any cyclic component. Choose a point fu(i,i) 

let f denote the lifting ofE^f. 
Then every x G S^{f^ I) has a unique lifting x 
w.r.t. the base point q fulfilling |^zy(i,i) — fu{i,i)\y,m < 
I. Furthermore, 

Mx) = Mx) for all xGS^{f,l), (36) 

O 

where Jq either denotes the functional from (14) or 
from (15) and Jq is its analogue in y by (35). 

Proof. Let us consider x G <S^(/, |). For {k,l) in the 
eight-neighborhood of {i,j), we have 


where 


eoo{x, f) := max dx,m{xij, (31) 

and the mapping 


dx-) ^k,l) — dx ,m ipi-tj 5 fu(i,j)) 


T dx,m{fu(i,j)^ fu(k,l)) 

Xdx ,m i,fy{k^l) •) ^k,l ) 


dl(f)<^. (37) 


u: Qq ^ assigns {i,j) G Qq 

a nearest neighbor in Q 


c(32) 


Here we measure the vicinity with respect to the dis¬ 
tance on the grid induced by taking shortest paths 
with respect to the eight-neighborhood. The x speci¬ 
fied this way are ‘near’ to the images / on and do 
not vary too much in Q. We also need an extension 
operator E extending a function / defined on to a 
function E{f) defined on Hq- A particularly simple ex¬ 
tension operator is the nearest neighbor operator Ej^, 
defined, for all (i, j), by 

E^fij = (33) 

with u as in (32). Then there is a constant inde¬ 
pendent of / but dependent on fl, such that 

dUEAf)) < BAmd^if). (34) 

We now consider the ‘lifted’ inpainting functionals. 
We first note that we may write the problem (14) in 
the form (15) modifying the distance term to be infi¬ 
nite if X 7 ^ / on . Then, on the universal covering 
space y of X, the inpainting functional Jq reads 

Jn{x) = FQc{x-,f)+ aTyi{x) 

^ ^ (35) 

+ ^TV2(a:)+7TVi.i(a:), 

where lifted image of /. The function- 

als F, TVi, TV 2 TVip are given by the correspond¬ 
ing functionals of Section 2.1 (there denoted without 
tilde) using y as the underlying vector space. We get 
the following relations: 


By assumption, we have dx,m{xuii,i)^ fuii,i)) < f • 
Therefore, every x G 5(/, |) has a unique lifting x by 
Lemma 3.6 w.r.t. to the base point q fulfilling \\xi^i — 

/mII <1- 

In order to show (36) we show equality for each 
of the involved summands. First, we consider TVi. 
By Lemma 3.6 we have dx,m{xij,Xk) = \xij — 
Xk\y,m^ k G {(i,j + l),(i-hl,:/Mi + l,j + l)}. Hence 
the definitions of TVi and TVi imply TVi(x) = 
TVi(:r). Concerning second order differences, we first 
consider the expressions D 2 {xi-ij, Xij, Xi^ij)- Simi¬ 
lar to (37), we have dx,m{xi-i,j^ < f which 

implies that the distance between any two members of 
the triple is smaller than Due to the properties of 
the lifting x this implies D{xi-i^i, Xij , Xixij]b 2 ) < tt. 
The same argument applies to D2{xij-i,Xij,Xijxi) 
which yields the equality for the TV 2 terms. A 
similar argument shows that TViq((r) = TViq(:r). 
Concerning the data term E{x] f) we define rij = 
dx-) fv(i,j)^ and r^j fi^(i,j)\y,'m‘ ^Ae 

show that rij and Xij agree for any (i,j) G 
By definition of S^{f,6), we have | for all 

(i,j) G flo- Furthermore, by the construction of / 
and X it holds hij = -h 27rkij, with kij G N and 
= 0 - We estimate |7i,j+i — ~ 

fv'{i,jA-l)\y,m~\d^i,j~fu{i,j)\y,m\ ^ 4 - ^i,j ki^xh 
then there exists k G ^\{0} such that 

|n,i+i - ri,j \ = Injxi - nj + 27rk\ > 27r - f > f. 

This is a contradiction and therefore kij = kij^i- 
Similarly we conclude kij = ki^ij- Hence, kij = 
^zy(i,i) = 0 foi* all (bi) ^ ^0 which implies rij = hij 
for all {i,j) G and completes the proof. □ 
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To formulate the next lemma we need the quantity 
di for functions defined on the complement of the in¬ 
painting region. We consider /: T, and define 

d^{f) in analogy to (30) by 


implies dp^^rn{x[jJij) = dx,m{xlj^ fid) all (i, j) G 
Qp. We extend the data / given on Qp to a grid func¬ 
tion g defined on CIq by setting gij = {Ej,f)ij = fv{i,j) 
for all {ip) G We get the estimate 


d^U)= E 




max dx 


{fi,j 5 fk,l): 


with Afij being again the eight-neighborhood of (iP). 
For the nearest neighbour extension operator Ej^ de¬ 
fined in (33) we have the following estimate: there is 
a constant independent of / but dependent on ft 
(and on f^o), such that 


d,iEM)) ■■= dfiE.if)) < C,iQ)d2{f). (38) 


As a further preparation, we need the following ob¬ 
servation. We consider the pure inpainting functional 
(14). We let u* be a minimizer of (14) for given data 
/. Then there is a constant B'pVt) > 1 which depends 
on 11, but not of /, such that 

dx,m{<vk 3 )<Bimd%{f). ( 39 ) 

For functions / with small values d^(/), this estimate 
follows from the boundedness of second differences by 
first differences. For the remaining /, the estimate 
(39) follows from the boundedness of the sphere as 
a set. For vector space valued data, the boundedness 
of second differences by first differences implies the 
estimate for arbitrary input. We note that we inten¬ 
tionally choose the symbol to avoid confusion when 
applying (39). 

Lemma 3.8. Let e > 0 and consider /: ft^ -d 
A” with d^{f) < 4 ,bi(q) ■ define p := 

max{(ai,..., 0^4, ydi, yd2,7}. and assume that p is so 
small that 


7 (/) < 


•20pcpn)Bpn)^ 


(40) 


where Cj^{fl) is given by (38). Then any minimizer x"" 
of the inpainting functional Jq given in (14) or (15) 
fulfills 

ecx)(^*,/) < 


Proof. For a minimizer x* of Jq given by (14) or 
by (15), we consider x' which we define as the clos¬ 
est point extension Ey of the restriction x*|qc to the 
non-inpainting region ft^. In view of the definition of 
600 in (31) we consider the estimate 

^ (^i^j 5 ) (41) 

We start with the second summand on the right hand 
side and notice that x'- j = xfj for (iP) G ft^ which 


MxP < Jo(^) =^TVi(^)+^TV2(^)+7TVi,i(^). 

(42) 

We further estimate the right hand side of (42): 
since the second order differences may be estimated 
by two times the first order differences, we get 
(3 TV2(^) < 2(max^ A)(l: 0? 0) TVi(^), and an anal¬ 

ogous inequality for 7TViq(^). Hence 

Jn{g) < 5max{(ai,... ,0^4, A,^2,7} ( 1 , 1 , 1 , 1 ) TVi(^). 


Next, we estimate each summand appearing in TVi(^) 
by the corresponding summand in dpg) to conclude 
that (1,1,1,1) TVi(^) < Adpg). We use (34) to get 

Ja{g) < 2t)pdi{g) < 2QpCu{U)d^{f). 

As a consequence we obtain for (f,j) G $7*^ 

dx,m{x'ij,fi,jf = dx,m{xlj, fi,jf < Jn{x*) 

< Mg) < 20pMCt)dkf) (43) 

< ( _g_i2 ^ /e',2 

By the definition of x', this implies that for all (ip) G 
LIq, we have 

,m{Xij 1 fu(i,j)) ^ 3* 

Next, we look at the first summand in (41). Since 
M M kj) e we may restrict to estimate 
dx,m{x'ij^x'l j) on ft. For this purpose, we consider the 
pure inpainting problem (14) for the inpainting region 
ft for data x*. Using (39) with data / = x* we have 
that 

dx ,m {Xij •) Xij ) 

< BUn) d^{x') 

< Bl{Q,){d^{f) + 2 max dxMkjJHhj))) 

<§• 

For the second inequality we used (43). Combining 
the estimates (44) and (45), we get 

dx,m{Xij^ fy(ij)) E dx,m{Xij^Xj^j) E dx ,m{Xij ^ fu(ij)) 

- I + I = ^ 

which implies that 

eoo{x\f) = ^ max {dx,m{xljju(i,j))} < e 
(i,j)GOo 

and this finishes the proof. □ 
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Lemma 3.9. The statement from Lemma 3.8 does 
also hold for data f : y and the inpainting fune- 

tionals Jq given by (35). 

Proof The statement is obtained following the lines 
of the proof of Lemma 3.8. □ 

Now we combine Lemma 3.7 and 3.8 to locate the 
minimizers of J and J. This part of the proof is rather 
similar to [ 6 ] which is the reason for streamlining it. 

Lemma 3.10. Consider e with 0 < 5 < | and let the 
funetion f: T fulfill both d^{f) < as 

well as d^{f) < and assume that the parame¬ 

ters (a, yd, 7 of Jq from (14) or (15) fulfill (40) w.r.t. 

6 . 

Then any minimizer x* of Jq lies inS^{f^ |). Fur¬ 
thermore, denote by f the unique lifting of f w.r.t. a 
base point q and fixed fv{ip) with 7rg(/^(i^i)) = fu(ip). 
Then eaeh minimizer of Jq defines a minimizer 
X* := 7rq{y'^) of Jq. Conversely, the uniquely defined 
lifting x* of a minimizer x* of Jq is a minimizer of 
Jo. 

Proof. By Lemma 3.8, any minimizer x* of the in¬ 
painting functional Jq fulfills eoo(^*,/) < Since 
5 < f, we get X* G S^{f, |). For the second state¬ 
ment we notice that the mappmg x i-G x is a bijec- 
tion from S^{f,6) to the set S^{f,6) := {y: Qq 
IV ij - Lii,j)\y,m < for every 6 with 0 < ^ < |. 
In particular, we may choose 6 = |. If G J’ is a 

minimizer of Jq, it lies in S^{f, |) by Lemma 3.9. 
By (36) and the minimizing property of y* we obtain 
for any x e 5“(/, f) that Jai^giy*)) = Jaiy*) < 
Jn{x) = Jn{x). As a consequence, 7rq{y'') is a mini¬ 
mizer of Jq on 5^(/, |). By Lemma 3.8 all the min¬ 
imizers of Jq are contained in S^{f, |) and hence 
7rq{y'^) is a minimizer of Jq. For the last state¬ 
ment let X* be_a minimizer of Jq. For its lifting x* 
and any y G S^{f,^), we get Jo(J*) = Jo(^*) < 
Jn{7Tq{y)) = Jn{y)‘ Thus, x* is a minimizer of Jq 
on S^{f, |). Since by Lemma 3.9 all minimizers of J 
he in <S^(/, |), the last assertion follows. □ 

After establishing relations between these function¬ 
als and their lifted versions, grid functions and data, 
we next formulate a convergence result for vector space 
data in A’. It is a reformulation of a convergence result 
which can be found in [ 2 ] for the more general class of 
Hadamard spaces or which can be derived from [ 8 ]. 

Theorem 3.11. Let J = YTi=i with eaeh Ji being 
a proper, elosed, eonvex funetional on and assume 
that J has a global minimizer. Assume further that 
there is L > 0 sueh that the iterates of the 

CP PA, ef. Algorithm 1, fulfill 


for all k G No. Then the sequence converges 

to a minimizer of J. In partieular, the iterates fulfill 

||^(^+^) _ II < 2XkL, (46) 

and, for all x G 

||^(^+i) — x||^ < 

-2Xk{J{x^^^)-J{x)) (47) 
+ 2A^Z/^c(c “h I). 

Next we locate the iterates of the CPPA for vector 
space data in A’ on a ball whose radius can be con¬ 
trolled. Since the data /: ^ A’ is not defined on 

the whole grid flo, we incorporate an extension oper¬ 
ator E, e.g. Ej^. An extension is needed as an initial¬ 
ization of the CPPA. We note that there is a positive 
number L' such that the iterates j produced 

by Algorithm I fulfill 

||^(/)-^(fe+i)| 7 <L', (48) 

This can be seen by taking mo,Mo G M as the min¬ 
imum and maximum of all components and pixels 
of E{f), respectively, and by letting G M 

the corresponding minima and maxima of the iterates 
x^^+c). Looking at the concrete form of the proxi¬ 
mal mappings, the only proximal mappings which can 
increase the maximal value or decrease the minimal 
value during the iteration are those of the second or¬ 
der differences. The possible increase is immediately 
decreased to initial niveau (of the macro-step) by the 
proximal mappings of data and first order terms if 
outside [C • mo,C • Mq] for sufficiently large C. We 
note that L' does not depend on the particular vector¬ 
valued / but only on Mq — mo. 

Lemma 3.12. Let f: y and a parameter se- 

quenee X = {A/c}/c of the CPPA with property (16) be 
given. Further let j be the sequenee produeed 

by Algorithm 1 for the inpainting funetionals Jq given 
by (35). Let x"": LIq ^ y be a minimizer of Jq. Then, 
for all /c G No and all I e {I,..., c}, we have 

- X*\\ < R (49) 

where 

R ■■= V||i;(/)-a:*||2 + 2||A||2L2c(c+l) + 2||A||oocL, 

where L = max(4, L') using L' from (48) and c denotes 
the number of inner iterations, i.e. c = 18 in ease 
of (14) and c = 19 in ease of (15), respeetively, and 
E is an operator extending y valued funetions defined 
on to fto used for initializing the algorithm. Here, 
7 Halloo = supj Ai. 
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Proof. Equation (47) in Theorem 3.11 tells us that 

-2Afe(J(a;('=)) - J(a;)] (50) 

T 2A^T^c(c T 1). 

We choose L" as the maximum of the Lipschitz con¬ 
stants for the terms Ji of Section_^3 origin^ng from 
the splitting of the terms TVi, TV2 and TVip. We 
note that these Lipschitz constants which are formed 
w.r.t. the Euclidean norm on y are all bounded by 4 
which, for the second order differences, is seen by esti¬ 
mating second order difference by sums of first order 
differences. This implies L" < 4. For the quadratic 
data term, we may differentiate Xij ~ 

Then we notice that the Xij are confined to an L' ball 
around E{f)ij (cf. (48)) which bounds the Lipschitz 
constant of the data term (where the metric on y is 
the Euclidean metric) by V in this case. Therefore, 
we can set L = max(I/',4). We apply (50) with a 
minimizer x = x* to obtain the estimate 


- a;* 


< \\x' 


< -a;*|| +2A|i:2c{c+l) 




i=o 


Using = E{f) yields < \\E{f) — 

-h 2||A|| ||2T^c(c + 1). Now we use (46) and the 
triangle inequality to estimate 


||^(^+c) — X*\\ 

< 2XkcL^\\x^^^ -x*|| 

< 2XkcL + ^/\\E{f) - x*P + 2\\XpL^c{c+ 1). 


This completes the proof. 


□ 


The following lemma shows that lifting commutes 
with applying the proximal mappings for the previous 
assumptions. 

Lemma 3.13. Let f: Qp T with d^{f) < 

and its lifting f of f w.r.t. a base point q as before. 
For eaeh summand Ji in the splitting Jq = Ji from 
Seetion 3.3 for both inpainting funetionals Jq defined 
via (14) and (15), their eorresponding funetionals Jq 
from (35), any x G S^{f, |), and its lifting x w.r.t. g, 
we have 


Prox^j, (x) = 7r,(prox,^j- (x)), (51) 


/ for indices {i,j) € . Since x € S^{f, |), it holds 

dx,m{xij,fij) < I for all (z,j) e . The com- 
ponents of the proximal mapping proxj^j^^ are given 
by Proposition 3.5 from which we conclude (51) for 
I = 19. The other proximal mappings of Ji,..., Jig, 
are given via proximal mappings of the first and sec¬ 
ond order cyclic differences from Theorem 3.3. We 
first consider cyclic components of the first order dif¬ 
ferences, i.e. summands involving Di. By the triangle 
inequality we have 


dx ,m {Xij ^Xij-^i) 

^ dx ^ fiy(i,j)) T dx ,m{fiy(i,j) f 

<^+B,(n) dS.(/)<f- 

Analogously we get dx,m{xi,j,Xi+ij) < By the 
explicit form of the proximal mapping given in The¬ 
orem 3.3 we obtain (51) for the J/, I = 1,...,8, 
which involve first order differences. Next we consider 
the second order differences D 2 with respect to the 
cyclic components. Let us exemplarily consider the 
vertical second order difference ^i,j+i)- 

Analogously as above, we see that the inequali¬ 
ties dx— li ^i,j^ “^5 dx 

and dx,m{xij-i^Xij) < ^ hold. Hence all the cyclic 
parts of the contributing values of x he in a common 
ball of radius 7r/2. Applying the proximal mapping in 
Theorem 3.3 the cyclic parts of the resulting points he 
in a common open ball of radius tt. 

An analogous statement holds true for the horizon¬ 
tal part. Hence the proximal mappings of these second 
differences agree with the cyclic version under iden¬ 
tification via TTq. This implies (51) for J9 ,...,Ji 4. 
It remains to deal with the cyclic components of the 
mixed second order differences Di^i{xij^ Xi-^ij^ Xij-^i^ 
Xi-^ij-^i). As above, we have for neighboring data 
items that the distance on the cyclic parts is smaller 
than ^. For all four contributing values of x we have 
that the pairwise distance on the cyclic parts is smaller 
than So their cyclic parts again lie in a boll of ra¬ 
dius smaller than tt and the proximal mappings agree 
under identification. This completes the proof. □ 

In the following main theorem we combine the pre¬ 
ceding lemmas to show that the output of the applied 
proximal mappings remains in 5^(/, |). This then 
allows for an iterated application of Lemma 3.13. 

Theorem 3.14. We ehoose the sequenee X = {A/c}/c 
fulfilling property (16) and 5 > 0 sueh that 


for all I G {!,..., 18} in ease of (14), and for all I G 
{!,..., 19} in ease of (15). 


y4£2 + 2||A||2L2c(c+l) + 2||A||c„cL < E 


Proof. The functional J19 from (29) appearing in the where c = IS or c = 19 and L = max(4, L') with 
splitting of (15) is based on the distance to the data L' as in (48). We eonsider data f: A with 
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both d^if) < and d^{f) < We as- 

sume further that the parameter veetors o^, 7 of the 

inpainting funetionals Jq given by either (14) or (15) 
satisfy (40) and that the initialization of the inpainting 
region E{f) is elose to the nearest neighbor extension 
Ey{f) in the sense that 

dx,m{E{f),Ei,{f)) = maxdx,miE{f)ij,E^{f)ij) < e. 

Then the sequenee generated by the CPPA 

given by Algorithm 1 eonverges to a global minimizer 
ofJ. 

Proof. Let / be the lifting of Ejy{f) with respect to a 
base point q not antipodal to fv{i,i) and fixed fv{i,i) 
with 'Kq{fv{ip)) = /zy(i,i)- Furthermore, let Jo denote 
the analogue of Jo for y valued data given by (35). 
Since dx^m{E{f),Ey{f)) < 5 < ^, the function E{f) 
is in 5^(/, ^) C iS^(/, |) (which is important for the 
application of Lemma 3.7 to the grid function E{f) 
later on.) From Lemma 3.8 we conclude that the min¬ 
imizer of Jo fulfills 11^*—/II <^< where / 
is defined at the beginning of this proof. By (49) we 
obtain 

R = + 2||A||2L2c(c+ 1) + 2||A||oocL 

< f2\\y* - ff + 2||/ - E(f)f + 2\\\\\n^c{c + 1) 

+ 2||A||oocL 

< y4£2 + 2||A||2L2c(c+l) + 2||A||c„cL < 

where E{f) denotes the lifting of E{f). By 
Lemma 3.12 the iterates of the CPPA for y- 

valued data fulfill 

Hence, combining these estimates yields — 

/II < 11^^^+-^ - ^*11 +lb* - /II < If. Estimating the 
discrete norm by the norm, we get from the pre¬ 
vious line that ||^^^^c) — /||oo < |f which means that 
all iterates stay within 5^(/, |). 

After these preparations we now consider the se¬ 
quence of the CPPA for A-valued data / 

with initialization E{f). We show that = 

TTqiy^^^c)). 

By definition and by Lemma 3.10, we know 

= E{f) = TTq{E{f)) = Trq{y^°'>). 

We continue a proof by induction and assume that 
^ = expg (^(^+ c )). By the local bijectivity of 
the lifting shown in Lemma 3.7, and since ^ G 


<S^(/, |), we conclude g <S^(/, |). 

By Lemma 3.13, we obtain 

= 7rg(prox^^j^7''+^))) 

By the same argument as above we have again 

^ 5(J^ 1). Finally, Theorem 3.11 tells us that 

x^^^ = 7Tq{y^^^) 7Tq{y*) as k ^ oo 

and by Lemma 3.10 we finally obtain that x* := 7rq{y'') 
is a global minimizer of J. 

□ 

4 Applications 

In this section we apply our algorithms to various im¬ 
age processing tasks. These are denoising in HSV 
space in Section 4.1, inpainting in both a noise-free 
as well as a noisy setting in Section 4.2 for both syn¬ 
thetic as well as real world data. Finally, we apply our 
algorithms for denoising frames in volumetric phase¬ 
valued data —in our case, frames of a 2D film— in 
Subsection4.3. Our approach is based on utilizing 
the neighbouring I frames to incorporate the temporal 
neighbourhood. The idea generalizes to arbitrary data 
spaces and volumes consisting of layers of 2D data. 

The algorithms were implemented in MATLAB The 
computations for the following examples were per¬ 
formed on a MacBook Pro with an Intel Core i5, 
2.6 Ghz and Mac OS 10.10.1. For all following ex¬ 
periments we set the sequence {Xk}k to = ff and 
set the number of iterations k to be 400 as a stopping 
criterion for Algorithm 1. 

4.1 Denoising nonlinear color space data 

As a first application, we consider denoising color 
space data. Various nonlinear color spaces have 
been considered in the literature; examples are luma 
plus chroma/chrominance based spaces such as YIQ, 
YUV and YDbDr and HSL type color space such as 
HSL, HSI or HSD. We here consider the HSV (hue- 
saturation-value) color space: the hue component is 
cyclic, the saturation and the value component are 
real-valued. 

We apply our algorithm for denoising combined cyc¬ 
lic and linear data to these x valued data. We 
compare the results with the usual approach using the 
linear RGB color space. For both spaces, we compare 
our approach on the product space with a model that 
denoises each channel separately. Finally, we com¬ 
pare the results of all these approaches under differ¬ 
ent noise models: we impose Gaussian noise on each 
component in RGB space. We impose Gaussian noise 
to saturation and value, and wrapped Gaussian noise. 
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u; oqjA^ c>2’ ^ 8’^ 32’ 

22.08 0 804. 23.43^ 0 867. 


Figure 5 . Denoising an image with independent wrapped Gaussian and Gaussian noise cr = | on each of the HSV chan¬ 
nels. The RGB-based approaches (b), (c) produce less colorful results than the HSV-based approaches (d),(e). In 
contrast to channel-wise denoising (d), the combined approach proposed in this paper (e) gets the object boundaries 
more properly and outperforms the other approaches in both PSNR ans SSIM. 


i.e. wrapping the normal distribution on the circle 
or in other words computing mod 27r after adding the 
noise, to the hue component of the HSV space. 

In Figure 5, a colorful drawing of a sailboat^ of 
size 512 X 512 pixel is obstructed by noise on all three 
channels of the HSV color space: for the hue, which 
is given on [0,1] we applied wrapped Gaussian noise 
(mod 1) with standard deviation cr = ^ and mean 0, 
for saturation and value, which are also given on the 
same range but are not cyclic, we applied Gaussian 
noise, also with cr = ^ and set all pixels exceeding 1 
to 1 and all deceeding 0 to 0. The resulting image 
is shown in Figure 5(a). We then apply four differ¬ 
ent first and second oder differences based approaches, 
where for each, the best result among the range of pa¬ 
rameter from a := ai =0^2 G ^Nq, i.e. the grid con¬ 
sisting of multiples of ^ including 0 as long as 7^ 0 
and the same range applies for ^ := /3i = /32 = 7. We 
further for this example 0^3 = 0^4 = 0, i.e. we focus on 
the purely anisotropic first order discrete total varia¬ 
tion and a second order isotropic version. We measure 
the quality using a peak signal to noise ratio (PSNR) 
and the structural similarity index (SSIM) on RGB 
color space and. We record the best result with re¬ 
spect to PSNR in 5 (b)-(e) and the original is shown 
in 5 (f) for comparison of the enlarged regions. 

First, we apply a real valued approach to each of 
the RGB channels separately. For a = P 
we obtain the best PSNR of 21.41, which is shown 


^Taken from the USC-SIPI Image Database, see http://sipi. 
use.edu/database/database.php?volume=misc&image=14 


in Figure 5(b). Applying a vector valued approach 
on RGB, i.e. setting m = 0, n = 3 in Algorithm 1, the 
best result obtains a PSNR of 21.46 for a = {3 = cf. 
Figure 5(c). While this outperforms the component 
wise denoising by taking combined color space edges 
into account, it does not reconstruct the colorfullness 
of the image. 

On the other hand, we apply Algorithm 1 to each 
of the channels of HSV, i.e., setting m = l,n = 0 for 
the first channel and taking the real valued case from 
above for the second and third. In order to keep the 
channels unsealed, the algorithm presented in this pa¬ 
per is rescaled to run w.r.t. to modi. The result is 
shown in Figure 5 (d) yielding a PSNR of 22.08. Fi¬ 
nally, applying a vector valued approach on HSV, i.e. 
setting m = 1 and n = 2 —again having the cyclic 
channel w.r.t. mod 1— yields an image shown in Fig¬ 
ure 5 (e), having a PSNR of 23.43, the best result of 
all four compared algorithms. Note that especially the 
colors are much better reconstructed than in the RGB 
based denoising approaches, which both suffer from re¬ 
duced saturation. This increase in quality can also be 
seen in the SSIM values denoted at the subfigure cap¬ 
tions in the bottom right. Interestingly the SSIM value 
of the componentwise approach on RGB is just a bit 
better than the vector valued approach. Furthermore, 
edges can be much better recognized in the vectorial 
approaches than in both channel-wise approaches, see 
especially the magnified region of the sail. 

Nevertheless there is are two tradeoffs. One con¬ 
cerns the overall amount of the parameters: Increas- 
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(a) RGB vectorial, (b) HSV vectorial, 
a = — B = — a = - B = — 

^ 32’ 16’ ^ 8’ 32’ 


21.42, 0.699. 23.43, 0.867. 


Figure 6. The best results on ^No of the parameters and 
the noise described in 6 but with respect to SSIM. 
While the best result w.r.t. SSIM coincides with the 
best PSNR value for the vectorial HSV approach in, 
see (b), the best vectorial RGB approach only increases 
slightly in SSIM, still suffering the same effects. 


ing the parameters smoothens the image and removes 
noise in larger constant regions like the sky but also 
destroys small details, like the leaves or the sail. The 
best PSNR is something in between: while the sky still 
resembles a little noise, the sail is only smoothened 
a little, but its main features are kept. The second 
tradeoff is between the first and second order terms. 
While a dominant value of a would keep edges in the 
image, it also introduces the well-known stair casing 
effect. This effect is reduced by the second order term, 
i.e. the parameter(s) which also smoothens edges. 
Due to the existence of both smooth regions and edges 
in natural images, a certain equilibrium has to be cho¬ 
sen. 


4.2 Inpainting nonlinear color space data 

Inpainting noise-free data. Here we consider the 
situation where some data items are missing, are lost 
or have been removed by a user. As example space 
we again consider the HSV color space. As in Section 
4.1, we compare the results with the usual approach 
using the linear RGB color space and with the HSV 
approach working component-wise. 

We consider a synthetic image in the HSV color 
space given by the function 

(atan2^, l-x^, l-|a; + y|), x,y & (52) 

where the first component is the arctangent function 
with two arguments. This extends a synthetic ex¬ 
ample used by the authors in [7]. The original image 
is shown in Figure 7 (f). The initial data is obtained 
by removing a disc with radius r = | as shown in 
Figure 7 (a). The goal is to “recover” the image in 
Figure 7 (f). For the inpainting we again apply Algo¬ 
rithm 1 using K = 800 iterations and performing a 
parameter search on g := gi = a2^ P •= Pi = P = 
2 = 7 G |No. Then, the real valued approaches in 
RGB color space, both channel-wise (b) and vectorial 
(c), do not reconstruct the original colors correctly. In 


contrast, the channel-wise approach (d) and the vec¬ 
torial approach (e) on the HSV space keep the colors 
and both produce very satisfactory results. The result 
of the vectorial approach in Figure 7 (e) is a nuance 
“sharper”, its reconstruction is a little bit better in 
PSNR at least. 

Inpainting and denoising data. In many situa¬ 
tions data is noisy and parts are lost or invalid. This 
results in an combined inpainting and denoising prob¬ 
lem for which we apply the proposed methods next. 
As in Section 4.1 before, we consider the HSV color 
space and compare the results with the usual approach 
using the linear RGB color space and with the HSV 
approach working component-wise. 

As a test scenario, we add wrapped Gaussian and 
Gaussian noise with a = | to the cyclic and non-cyclic 
components, respectively, similar to Section 4.1. Fur¬ 
thermore, we remove a ball in the center as also done 
in Section 4.2; see Figure 8. This initial is shown in 
Figure 8(a). Again, we would like to get back the 
image shown in Figure 8(f), which is data defined 
in HSV space and therefore also has features e.g. in 
the hue. The parameters for the algorithm are ob¬ 
tained using the same setup for Algorithm 1 as in Sec¬ 
tion 4.1. Then, the real valued approaches in RGB 
color space, both channel-wise (b) and vectorial (c), 
do not reconstruct the significant features. In con¬ 
trast, the channel-wise approach (d) and the vectorial 
approach (e) on the HSV space reconstruct all signifi¬ 
cant features and produce almost perfect results. The 
result of the vectorial approach in Figures 8 (e) is a 
nuance “sharper” yielding a slightly higher PSNR. 


Inpainting real world data. Looking at real world 
data, we compare an inpainting of the image “beach” 
In Figure 9 ((a)) some areas of the image were lost. We 
compare the vectorial approach on RGB color space in 
Fig. 9 ((b)) and HSV in Fig. 9 ((c)) for the same set 
of parameters a = P = ^: the HSV color model 
is able to also incorporate color changes as shown in 
the magnification, where the HSV model is able to re¬ 
construct the violet part of the clouds. We note that 
we employed the inpainting algorithm as described in 
Algorithm 1 also for the diagonal differences. Even it¬ 
erating alternatingly the anisotropic differences yields 
preference of the diagonal directions, which is even 
enhanced by the diagonal differences. Therefore both 
models tend to create wavy structures, when the sur¬ 
rounding is not isotropic as in the last examples. 

Despite the most bottom left inpainting area, the 
HSV model performs better than the RGB model with 
respect to PSNR, where the HSV model obtains 30.612 

^Available at http://pixabay.com/de/strand-lagune- 
sonnenuntergang-164288/ and published under public 
domain. 
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(c) HSV channels, (e) HSV vectorial, 

c^ = i/3 = i, « = = 

PSNR: 27.37. PSNR: 27.68. 


Figure 7. Reconstruction of a synthetic image with the black inner circle missing in a noiseless setup. The RGB-based 
reconstructions (a),(b) both yield a degrading of the colors, while the HSV-based approaches (c),(d) reconstruct 
the colors. The proposed approach (d) based on a x model yields the best PSNR. 


and the RGB model is slightly behind with 30.548. 
The same holds for the SSIM, though there the dif¬ 
ference is quite small, because even the SSIM of the 
lost version —setting lost pixel to white as in Fig¬ 
ure 9 ((a))— yields an SSIM value of 0.9202. We 
obtained an SSIM value of 0.98866 for the complete 
image with the HSV based inpainting and 0.98862 for 
the RGB based inpainting. 

4.3 Denoising sections in volumetric cyclic 
data 

Finally, we apply our algorithms for denoising frames 
in volumetric data. Examples of volumetric data are 
frames of a 2D film or a stack of slices, each slice being 
a 2D image as appearing, e.g., in computed tomogra¬ 
phy. We want to denoise such slices incorporating the 
temporal/spatial information stemming from the third 
dimension, which is, e.g., the temporal neighborhood 
information in a film. 

To be more precise, we consider volumetric data 
Ikihj)^ where is a vector valued pixel at the 

pixel location i^j in the kth frame/slice. The setup 
is rather general, and we can assume being 

data from some rather general space A4 - say a man¬ 
ifold; here we exemplarily consider = S^. We take 
a look at the I neighboring (left and right) frames 
Ih^h = /c — /c,..., /c + /, around a center frame 

/ = J/c at position k. With the position k, we now 
associate bivariate data living in being the vec¬ 

tor of data points at the same position (i,j) in the 
neighboring frames. To be precise, we consider the 
bivariate data J/^, with Jk{hj) ^ given by 


JkiiJ) = (4-ft(*,j'')) _ • We apply our algorithms 

to the derived data Jk and compare the result to the 
usual denoising of the single frame J/^. 

The video underlying Figure 10 is constructed as 
follows. As basis, we use the image given by the first 
component of (52) on [—|, outside the disc of ra¬ 
dius we add J which is the same as rotating the 
input (i,jf) of each pixel in Dq by the same amount 
clockwise. The video consists of 13 frames rotating the 
disc clockwise by tt and the outer region by tt coun¬ 
terclockwise, i.e. from before the center frame to 

at the end of the sequence for the inner and with 
changed signs for the outer redgion. This corresponds 
to a rotation of ^ per frame and region. We finally 
sample each of these frames with 256 x 256 pixel on 
r_i 112 
L 2 ’ 2j • 

In our example. Figure 10, we show in (d) the sev¬ 
enth frame of the constructed video from the previous 
the paragraph. On each frame, we impose wrapped 
Gaussian noise with standard deviation cr = |, see 
Figure 10 (a) for frame 7 of the video. In (b), we per¬ 
form denoising just on the frame k = 7. Performing 
a first and second order denoising yields staircasing 
and/or reduction of the sharp edge at the disc border. 
Choosing the parameters a = ai = a 2 ^ P = Pi = P 2 
from the set G and setting K = 400 as maximal 
number of iterations of the CPPA, we obtain the opti¬ 
mal value (a = ^, = 0. This indicates that staircas- 

ing still resembles a better result than unsharpening 
the edge, which would be the effect of including sec¬ 
ond order differences. In (c), we perform a combined 
vectorial denoising in (§^)^^ as proposed for a total of 
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(a) Synth, image 
with noise, cr = ^, 
lost area marked black 



(d) HSV channels, 

»=hP=l 

PSNR: 26.20. 


(e) HSV vectorial, 

a= !./?= i 

PSNR: 26.38. 


Figure 8. Reconstruction of a synthetic image with the black inner circle indicating missing data and the measurement 
itself being noisy. The RGB based reconstructions (b),(c) miss the main smooth features. The HSV model based 
reconstructions (d),(e) reconstructs these features yielding a satisfactory result. The proposed vectorial approach 
(e) yields the best PSNR. 



(a) Lost area (white). (b) RGB model. (c) HSV model. (d) Original image. 

Figure 9. Inpainting the real world image “beach”: The HSV based model better reconstructs the clouds than an 
RGB based approach. 


13 frames (/ = 6, /c = 7) and show the central seventh 
frame. At the cost of being computationally more ex¬ 
pensive due to the increased data set, this approach 
outperforms the first, single frame based approach. 

Discussion and Future Research 

In the following we discuss the relation of the present 
work to the author’s previous work [6, 7]. These works 
consider the -valued situation, whereas we here con¬ 
sider the product space x These product 

spaces are practically relevant: they appear, e.g., in 
the context of color space and, more abstractly, when¬ 
ever considering polar coordinates (magnitude and 
phase). In contrast to general manifolds the spaces 
considered here still allow for a solution to the con¬ 


sidered minimization problem, i.e. the involved prox¬ 
imal mappings, to be given in explicit form. We are 
convinced that this is no longer possible for general 
manifolds. The product spaces further imply changes 
in both the algorithms and the convergence analysis. 
Although this paper as well as [6, 7] (and also the 
previous work [85]) use a CPPA, the derived proximal 
mappings are different. In particular the mappings de¬ 
rived here do not arise as component-wise application 
of the one dimensional situation. This is especially 
not the case for the setting, where several data items 
are fixed by constraint. In this context, we again men¬ 
tion that the natural second differences on (§^)^ x 
which we defined in this paper couple the components 
of the range space, which also affects the convergence 
analysis The employed methods are based on the ones 
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(a) center video frame, (b) Single frame, (c) Vectorial frames, (d) Original Frame. 

Gaussian noise, a = a = PSNR: 29.69. ck = ^5 PSNR: 35.01. 


Figure 10. Denoising S^-valued video frames. The proposed approach for denoising S^-valued video frames (c) yields 
signihcantly better results than frame-wise denoising (b). 


used in [6] (which themselves are based on [82]), but 
the concrete analysis involves an additional degree of 
complexity stemming from the mentioned coupling. 

This work also contains material which is even new 
for -valued data. First, in contrast to pure de¬ 
noising, additional/different proximal mappings are 
needed when dealing with the inpainting problem due 
to the additional constraints. In [7], we simply used 
the proximal mappings computed in [6] and applied 
projections to ensure the constraints. Here, we com¬ 
pute the proximal mappings of the constrained prob¬ 
lems explicitly. This is more natural and also advanta¬ 
geous in the analysis. Concerning the analysis, we here 
include an inpainting setup. Even for the setting, 
this has not yet been done. In the conference proceed¬ 
ing [7] we do not provide a convergence analysis. The 
paper [6] which provides an analysis considers func¬ 
tionals for denoising in the setting only. 

A topic of future research are algorithms for higher 
order TV-type functionals for data living in more gen¬ 
eral manifolds. 
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